# ------------------------------------------------------------------------------
# Compares MACE rates and yield rates by
# Author: Cassidy Shubatt <cshubatt@gmail.com>
# To run: bash 04_mace_vs_yield.sh
# ------------------------------------------------------------------------------

# Libraries --------------------------------------------------------------------
library(here)
library(yaml)
library(data.table)
library(tidyverse)
library(glue)
library(testit) # other()
library(ggplot2)

u <- modules::use(here::here("lib", "util.R"))
a <- modules::use(here::here("lib", "aesthetics.R"))
a$get_font("Optima", here("lib", "optima.ttf"))

overnight_lab <- ""
temp <- here("code", "04_cost_effectiveness", "temp")

# Load Data --------------------------------------------------------------------
message("Loading data...")
paths <- read_yaml(here::here("lib", "filepaths.yml"))
demographics <- readRDS(paths$analysis$demographics) %>%
  select(ed_enc_id, death_030_day)
test_cohort <- readRDS(glue(paths$analysis$test_cohort)) %>%
  filter(!exclude) %>%
  u$safe_left_join(demographics)

# Ventile Rates ----------------------------------------------------------------
message("Getting outcome rates by ventile...")
outcome_rates <- test_cohort %>%
  group_by(tile_stent_or_cabg_020_tested) %>%
  summarize(
    num_unt = sum(!test_010_day),
    test_rate = mean(test_010_day),
    yield_rate = sum(stent_or_cabg_010_day) / sum(test_010_day),
    adverse_total = sum((macetrop_030_pos | death_030_day) & !test_010_day),
    adverse_rate = adverse_total / sum(!test_010_day)
  ) %>%
  ungroup() %>%
  select(-adverse_total)

# Plot -------------------------------------------------------------------------
message("Plotting mace rate vs. yield...")
# this number comes from 06_cost_effectiveness/02_cost_effectiveness_stats
risk_threshold_150K <- readRDS(paths$analysis$costeff_risk_cutoff)
fit <- lm(adverse_rate ~ yield_rate, outcome_rates)
corr_mace_threshold <- predict(
  fit,
  newdata = list(yield_rate = risk_threshold_150K)
) %>%
  .[[1]] %>%
  signif(3)

gg <- ggplot(
  outcome_rates, aes(
    x = yield_rate, y = adverse_rate, color = test_rate, size = num_unt
  )
) +
  geom_point() +
  geom_smooth(method = "lm", color = a$disc_palette[1], fill = a$disc_palette[1], alpha = 0.1) +
  geom_vline(xintercept = risk_threshold_150K, linetype = "dashed", color = a$disc_palette[1]) +
  geom_hline(yintercept = corr_mace_threshold, linetype = "dashed", color = a$disc_palette[1]) +
  annotate(
    "label",
    x = risk_threshold_150K, y = 0.1, family = "Optima",
    label = "150K cost-effectiveness", size = 14
  ) +
  annotate(
    "label",
    x = 0.2, y = corr_mace_threshold, family = "Optima",
    label = glue("MACE threshold: {corr_mace_threshold}"), size = 14
  ) +
  labs(x = "Yield of Testing", y = "MACE Rate", color = "Test Rate") +
  theme_bw() +
  scale_color_gradientn(colors = a$ordered_palette %>% rev()) +
  theme(
    legend.position = "bottom",
    text = element_text(family = "Optima", size = 40)
  ) +
  guides(size = FALSE)

# Save
message("Saving...")
ggsave(file.path(temp, "mace_vs_yield.png"), width = 10, height = 7, unit = "in")

message("Done.")
